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1 Introduction 

A solid transports energy. Besides the mobile electrons, one important mechanism 
for energy transport are the vibrations of the crystal lattice. There is no difficulty 
in writing down the appropriate lattice dynamics. To extract from it the thermal 
conductivity remains a fairly untractable problem. The most successful approach 
exploits that even rather close to the melting temperature the typical deviations 
of the crystal atoms from their equilibrium position are small as compared to the 
lattice constant. This observation then leads to the phonon kinetic equation, which 
goes back to the seminal paper by Peierls pQ. (For electron transport a correspond- 
ing idea was put forward by Nordheim [2].) Phonon kinetic theory flourished in 
the 50ies, an excellent account of the 1960 status being the book by Ziman [3]. Of 
course, transport of heat and thermal conductivity remain an important experimen- 
tal research area, in particular since novel materials become available and since more 
extreme properties are in demand. On the other hand, if the very recent collection 
of articles by Tritt [4j is taken to be representative, it is obvious that after 1960 
hardly any new elements have been added to the theory. The real innovation are 
fast and efficient molecular dynamics algorithms. The currently available techniques 
allow the simulation of 6 x 6 x 6 periodized lattices with two atoms per unit cell [5] . 

According to the Green-Kubo formula the thermal conductivity is determined 
through the time-integral over the energy current correlation in thermal equilib- 
rium. In my contribution I will explain its structure for weakly anharmonic lattices. 
While I do not add anything novel in substance, I believe that, with the post 1960 
insights gained from the kinetic theory of rarified gases, the story can be presented 
more concisely and systematically than done usually. As a bonus, the mathematical 
physics issues left unresolved will become more sharply in focus. 

2 Anharmonic lattice dynamics 

Physically, one starts from a given crystal structure, which means to specify the 
lattice and the number of atoms per unit cell. The interaction potential is expanded 
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in the displacements away from the equilibrium positions. Then the first order term 
vanishes, because one expands at a stationary point. The second order term is the 
harmonic approximation and higher order terms are regarded as small corrections. 
It is argued that for real crystals mostly the third order term suffices unless there 
are special symmetries which make it vanish and requires to go to fourth order. In 
this article, the focus will be on the analysis of the linearized Boltzmann equation 
and its relation to the energy current correlations. For this purpose we take the 
liberty to employ a single band model for the anharmonic lattice dynamics. There 
is no difficulty, in principle, to add on extra features so to make the model more 
realistic. 

We assume a simple hypercubic lattice Z d with a single atom per unit cell. 
Physically d — 3, but we keep the general dimension d because of recent interest 
in chains, for which d — 1. A single band model corresponds to scalar atomic 
displacements. 

Fourier transform will be convenient. Let T d = [— |, be the first Brioullin 
zone of the dual lattice. For / : Z d — > R its Fourier transform, /, is defined by 

f(k) = J2^ 2nk - x f x . (2.1) 

Here k G T d and f(k) extends periodically to a function on M. d . The inverse Fourier 
transform is given by 

f x = f dke^ k x f(k) . (2.2) 

JT d 

For x G Z d the deviation away from x is denoted by q x G K. The corresponding 
momentum is denoted by p x G K. We choose units such that the atomic mass equals 
one. The harmonic approximation to the interaction potential reads 

f/harm(g) = ^ Yl a ( X ~ ■ (2.3) 

The elastic constants a(x) satisfy 

a(x) = a(-x) , \a(x)\ < 7 e" 7l|:c| (2.4) 

for suitable constants 70,71 > 0. Mechanical stability requires 

a(k)>0. (2.5) 

In addition, because of the invariance of the interaction between crystals atoms 
under the translation q x -w q x + a, one imposes 

^2a{x) = 0, i.e. 3(0) = 0. (2.6) 
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For an optical band, because of the internal structure of the unit cell, the condi- 
tion (12.61) is not satisfied, which in the framework of our model can be interpreted 
as adding to the physical harmonic interaction satisfying (12.61) a harmonic on-site 
potential of the form 

t/site(g) = ^o 2 E^ ( 2 - ? ) 

2gZ d 

The harmonic lattice dynamics is governed by the hamiltonian 

= ~ Yl (p£ + + \ Yl a ( x ~ y } qxq y ( 2 - 8 ) 



and has plane wave solutions with dispersion relation 

uj{k) = {ujl + a{k)) 1 ' 2 . (2.9) 

Clearly, oj{k) = u(—k) and u(k) > uj q > 0. We concatenate q x and p x into a single 
complex-valued field a(k) as 

a(k) = -L(y^(kjq(k) + i— L=p(k)) (2-10) 
v2 \/u)(k) 

with the inverse 
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q(k) = —-^==(a(k) + a(-ky) , p(k) = y^{k) ( - a(k) + a(-k)*) . (2.11) 
The a-field evolves as 

d 

— a(k, t) = -iuj(k)a{k, t) . (2.12) 

In nature lattice vibrations are quantized. In our model this is easily implemented 
by promoting a(k)* and a{k) to creation and annihilation operators of a scalar Bose 
field. a(k)* is the operator adjoint to a(k) and the a(/c)'s satisfy the canonical 
commutation relations 

[a(k), a{k')*} = 5{k - k') , [a(k), a(k')] = . (2.13) 

The Heisenberg evolution for the a-field is still governed by (12.121) . 

Continuing the expansion scheme we add to Hq the next order terms. The 
simplest one would be a cubic on-site potential as 



3 

which in terms of the a-field reads 



^3 = ~5>*, (2.14) 



3 

V 3 =l I dk 1 dk 2 dh5{k 1 + k 2 + h)f\{2uj(k j ))- 1 / 2 {a{k j ) + a{-k ] )*) . (2.15) 

3 Jjad " 
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Correspondingly, at fourth order, 

V, = \ J2 4 (2-16) 



4 

x&L d 



which in terms of the a-field becomes 



i r 4 

V A = - - & 1 dMMfc4^i+^+fc3+^) f](2w(y )" 1/2 (a(y+a(-y *) . (2.17) 
4 Jjid J=i 

Hha, + AV3 is not bounded from below. This can be remedied by adding A 2 V4, for 
example, which would then not contribute on the kinetic scale. 

If the potential depends only on the displacement differences, then the lowest 
order nonlinearity is 

x,y& d 

with a 3 (— x) = — a 3 (x) and | ck 3 | exponentially bounded. Switching to the a-field 
V3di becomes 



Vm\ = 7. [ dk 1 dk 2 dk 3 S(ki + k 2 + k 3 

6 Jj4d 



x E a ^ H(Mk 3 ))- 1/2 {e^ x - 1) (afa) + a(-%)*) . (2.19) 

x& d J=l 

In the kinetic limit the square of the vertex function determines the collision rate. 
Thus, from the collision rate 

3 

n^^-))- 1 (2.20) 
i=i 

for the on-site V3 the collision rate for V^u is obtained by the replacement 

3 3 

n(2 W (y)- 1 |^a 3 WlI(e !2 ^-l)| 2 . (2.21) 

i = l x£Z d j = i 

Because of such a simple substitution rule we continue to work with V3. The corre- 
sponding rule also applies to the switch from V4 to V^i- 

From other areas of mathematical physics one is accustomed to have a given 
starting hamiltonian. In our context this means to specify the elastic constants a(x), 
a 3 (x), a.i(x). For real crystals their determination requires a lot of experimental 
(and also theoretical) efforts, as discussed in [4], see also 0, 0, [8] for a modeling of 
aluminium and silicon. It would be thus of importance to have a stability result 
available, which ensures that certain qualitative properties do not depend so much 
on the specific choice of elastic constants. 
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3 Energy current correlations 



Let us consider the Hamiltonian 

H = iJ ha + AF 3 + A 2 V 4 . (3.1) 

The total energy current correlation function is computed in thermal equilibrium 
at inverse temperature (3. It is denoted by C\(t) and will be defined below. Since 
A C 1, the plan is to compute C\(t) in the limit of A — > 0. The phonons then hardly 
interact and C\(t) decays slowly on the time scale A -2 . Thus one expects that the 
limit 

\imC x (\~ 2 t) = C kin (t) (3.2) 

exists and is determined by the phonon Boltzmann equation linearized at equilib- 
rium. 

Let us first find the local energy current. Since H is not local, there is some 
arbitrariness involved in defining the local energy. One conventional choice for the 
energy at site x is to set 

H x = \vl + \ u l<fc + ^^a(x- y)q x q y + ^Xq 3 x + ^\ 2 q* . (3.3) 

y&A 

In the Heisenberg picture H x becomes time-dependent. Writing dH x (t)/dt as a 
divergence, the energy current can be identified with 

J x = \Yl y a (y)(~ ( lxPx+y + Qx+yPx) (3.4) 
y£Z d 

which happens to be independent of A. To verify (13.41) . one chooses a large box A 
with faces <9A. The energy inside A is 

H h (t)=Y J H x (t) (3.5) 

and it satisfies 

j t H A (t) = J S) + 0(dA) , (3.6) 

xt=A 

where n x is the outward normal to A at x G dA. The errors come from the corners 
of A and from the possibly infinite range of a. 

With this input the total energy current correlation is defined by 

t ■ C x (t)£ =J2(( £ - J °W)( £ ■ J *(°))>/? > ( 3 - 7 ) 

x&Z d 

where £ G M d , J x G M d , " • " is the scalar product in R d , and C\(t) is a d x d matrix. 
(•)p refers to the thermal average with respect to Z~ x e~P u . By time-stationary and 
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time-reversal, C\(t) is symmetric and it suffices to consider its numerical range. At 
the expense of an error of order A, we may replace in e~@ H the full Hamiltonian H 
by the harmonic approximation 



H ha = / dku(k)a(k)*a(k) . (3i 



For the total current one finds 

J = V J x = — [ dk(Vu(k))uj(k)a{k)*a(k) , (3.9) 

where it is used that VS = 2oj(Vuj). Since [i?ha> J] = 0, ^2 x&Z d £ • J x can be lifted 
to the exponent. Thus we define the new average {-)p, T with respect to the state 
Z~ l exp [-(3H h& + r£-J]. Then 

£ ■ C x (t)£ = lim -(£ ■ J (t))p, r + 0(A) . (3.10) 

The anharmonicity now resides only in the dynamics. 

The limit A — > on the right hand side in (I3.10p is discussed in [9]. The initial 
state is spatially homogeneous and determines the Wigner function Wp iT (k) through 

(a(k')*a(k))p jT = S(k - k')W{,, T (k) (3.H) 

with i 

Wp )T (k) = (exp [Pu(k) - t(£ ■ Vw(k))w(kj\ - l) . (3.12) 

On the kinetic time scale, X~ 2 t, the Wigner function Wp >T evolves to W T (t) which is 
determined as the solution of the spatially homogeneous Boltzmann equation. For 
our model, i.e. for the anharmonic on-site potential V3, it reads 

^W{t)=C{W(t)) (3.13) 

with the collision operator 
C(W)i = - dk 2 dk 3 (uj 1 u 2 u z Y 1 

I Jf2d 

x {25^ + uj 2 - uj 3 )5(h + k 2 - k 3 ) (W!W 2 W 3 - W X W 2 W 3 ) 

+ b(uo x -U2- oJ 3 )6(h -k 2 - k 3 ) {W{W 2 W 3 - W X W 2 W Z ) } • (3.14) 

Here we use the shorthands Wj = W(kj), ojj = to(kj), j = 1,2,3, and W{k) = 
1 + W(k). Using (13.91) . the average in (13.101) becomes then 



(£.J (\- 2 t))p jT = ^- [ dk(£-Vu{k))u{k)W T {k,t) + 0{X). (3.15) 
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The next task is to take the limit r -> in fl3~T0l . One has Wp >0 (k) = (e^ k) - l)" 1 
which is a stationary solution of (13.131) . Thus the limit r — > amounts to linearize 
(I3.14p at the equilibrium Wigner function 

W^Jfe) = (e^W - I)" 1 , (3.16) 

to say 

W P;T = W P + rWpWp{l ■ Vuj)uj + C(r 2 ) . (3.17) 

Note that J dk(Vu)uW/3 = 0. As suggested by (I3.17p . with a significance which 
will become more convincing in the context of the Gaussian fluctuation theory, see 
Section [5], the natural linearization of C is 

C(W P + 5WpWpf) = -SLf + 0(8 2 ) . (3.18) 

From (I3.14p one deduces 




dk 2 dk 3 (uiu 2 uj 3 ) 



x (28{u x + u 2 - u 3 )5(h + k 2 - h)W pi Wp 2 W^{f x +f 2 - ; 3 ) 

+8(u! -U2- u 3 )8{h -k 2 - k 3 )W 01 W 02 W P s(fi -h~ h)) • (3.19) 

Properties of L will be discussed in the subsequent section. 

Let A be the linear operator obtained from flat linearization as 

C{Wp + 5f) = 5Af + 0{5 2 ) . (3.20) 

Clearly A(WpWpf) = -Lf. Combining fl3TT0l and fl3~7L5D we finally conclude 

\im£- C x (\- 2 t)e = £ ■ C ki Jt)£ (3.21) 

with 

£ . C kin (t)£ = {{2-k)-\£ ■ , e-^WpWp^ir)- 1 ^ ■ Voj)uj) , (3.22) 

where (•, •) is the inner product in L 2 (T d , dk). 

For future use it will be convenient to write Ckm(^) i n a more symmetric form. 
Expanding the exponential one notes that 

£ ■ C kin (ty = ■ Vuo)uo{WpWp) l/2 , 

exp [- (T^T^)- (1/2 >Z.(Vl^ 3 W /3 )- C:L/2) (W /3 W /3 ) 1/2 (2 7 r)- 1 (^ - Vc^)c^> . (3.23) 

As will be shown, L = L*, i.e. L is a symmetric operator in L 2 (T d ,dk). Therefore 
Ckin(^) is a positive symmetric d x d matrix. 
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In the kinetic limit the thermal conductivity is given through 

£n kiQ -£ = (3 2 / dt£-C km (t)£ 
Jo 

= (2tt)- 2 (3 2 ((£ ■ Vu)uWpWp , L-\£ ■ Vu)uWpWp) . (3.24) 

Reversing our argument, and assuming uniformity in t for the limit A — > 0, one 
infers that the true thermal conductivity, k(X), of the anharmonic model behaves as 

k(A) \~ 2 K kiQ (3.25) 

for small A. 

In the classical limit [a(k) , a(k')*] = 0, i.e., W = W. In the definition of L one 
has thus to replace 

W P , Wp by Wf{k) = . (3.26) 

We presented the argument for a cubic on-site potential. But, clearly, the result 
holds also for other small anharmonicities. Only the collision operator, and its 
linearization L, would have to be modified. 



4 The linearized collision operator 

If one accepts the argument leading to (13.211) . the remaining task is to study the 
spectral properties of the linearized collision operator, from which the time decay 
of Ckin(t) can be infered. While this looks like a conventional mathematical physics 
problem, the difficulty comes from the energy- momentum constraint. Only in a 
few special cases there is an explicit solution. Otherwise one has to work with the 
implicit definition. In fact, there can be no solution at all, in which case L = 0, or 
several solutions, in which case one has to sum over all collision branches, 
(i) quadratic form. For three phonon processes, on-site potential V3, the quadratic 
form of the linearized collision operator L = L3 is given by 

(g, L ? J) = - / dk 1 dk 2 dk 3 (uj 1 uj2^3y 1 5(uj 1 + uj 2 - ^3)^(^1 + k 2 - h) 

l Jj3d 

xWp 1 Wp 2 Wp 3 (g 1 + g 2 - g 3 )(f 1 + f 2 - f 3 ) . (4.1) 
Correspondingly for the on-site potential V4 one has 

3tt 3 
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J dkidk 2 dk 3 dk^{uoiUJ 2 uj 3 uJi) 1 5{uo\ + co 2 — uo 3 — u&)6(ki + k 2 — k 3 — k^ 



•x.Wp X Wp 2 Wp 3 Wp±{g x + g 2 -g 3 - # 4 )(/i + / 2 - / 3 - U) 
+ — / dkidk-zdkz.dk^ijJiuj-iUJ^i)' 1 5{u}\ + uj 2 + uj 3 — u}^)8{ki + k 2 + k 3 — k^j 

p X W p 2 W p 3 W p±{g x + g 2 + g 3 - ^ 4 )(/i + h + h ~ h) 
(f,L ip f) + (f,L 4t f). (4.2) 
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L 4p corresponds to the collision of a pair of phonons and L 4t to a merger of three 
phonons into a single one, and its time reversal. The quadratic forms for L 3 , L 4p , L 4t 
(for notational simplicity from now on commonly denoted by L) are somewhat for- 
mal. Firstly, if cu(0) = and u(k) > otherwise, the smooth functions g, f have to 
vanish at k — 0. More seriously, the proper definition of the 5-function requires to 
study more carefully the solutions to the energy constraint 

u{k x ) +co(k 2 ) = uj{h + k 2 ) , (4.3) 

say in the case of L 3 . For the purpose of our exposition, let us simply assume 
that the quadratic form defines L as a self-adjoint operator. Clearly, L > since 
(/, Lf) > 0. As (£ ■ Vuj)u) is bounded, one has 

I ■ C kin (t)£ < (2tt)- 2 ((£ ■ Vuj)u, {£ ■ Vw)w) . (4.4) 

(ii) zero subspace. To establish that lim^oo C^it) =0, (£■ Vw)w has to be orthog- 
onal to the zero subspace of L. There seems to be no cheap argument and one has 
to study the solutions to 

Lf = 0. (4.5) 

From ( 14. ip . ( 14. 2 p it follows that / has to be a collisional invariant, see P,[in] for the 
definition. Considering only the first summand of ( 14.21) . there is a general argument 
[lOj, that the solutions to (/, L 4p /) = are spanned by 1, uj. Note that (1, {£ ■ 
Vw)w) = 0, (u, (£ ■ Vuj)u) = 0. The constant function results from phonon number 
conservation in a pair collision. This conservation law will be destroyed by adding a 
little bit of either three-phonon, L3, or the second term of the four-phonon processes, 
L 4t . The zero subspace is then one-dimensional and spanned by uj only. For L3 of 
(14. ip . the classification of the collisional invariants is an open problem. 

(iii) spectral gap. If L has a spectral gap, the energy current correlation decays 
exponentially. If in addition (£ ■ ^7u)u is orthogonal to the zero subspace of L, 
then the conductivity, as the time-integral over Ckm(t), is finite (and non-zero). In 
particular ( I3.25P holds. 

L is a sum of a multiplication operator and an integral operator, 

L = V + I, Vf(k)=V(k)f(k), If(k)= f dk'l(k,k')f(k'), (4.6) 
where, say in the case of L 3 , 

v(k) = -w p {k)uo{k)- 1 [ dkifofaMk + h)y l 

2 Jjd 

x (25(uj(k) + w(ifei) - oo(k + h^Wpi^Wpik + h) 
+ 5{uo{k) - u(h) - u(k + k^Wpik^Wgik + kSj . (4.7) 
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The integral kernel I(k, k') is implicitly defined. It has no definite sign and tends to 
be divergent on lower-dimensional submanifolds of T d x T d . It would be useful to 
know under what conditions the integral operator I is compact. 

In the very common relaxation time approximation, / is simply dropped and one 
sets in approximation 

£ . C kin (t)£ = (2tt)- 2 ((£ ■ Vuj)oj, e- w/T (£ ■ Vw)w) (4.8) 

with the relaxation time 

r(k) = Wp(k)W p (k)V(ky l , (4.9) 

see (3.23). 

(iv) FPU chains. The Fermi-Pasta-Ulam chain is the special case d = 1 with near- 
est neighbor coupling and no quantization. For a harmonic on-site potential the 

1/2 

dispersion relation is u (k) = (wg + 1 - cos(27rJfe)) ' , k E T. Although d = 1, 
the conservation laws of energy and momentum allow for non-degenerate pair col- 
lision and L^ p ^ 0, while L 4t = [TTJ [121 13] • There are fairly explicit formuli for 
the potential V and the integral kernel / [H]. For ujq > and a quartic on-site 
potential V4, the linearized collision operator has a gap and the zero subspace is 
two-dimensional. The gap seems to close as ujq — > 0. On the basis of numerical 
simulations, the conductivity should be finite even for ujq = [12J. The FPU-/3 
chain has the nonlinearity V^. V and I has been computed by Pereverzev [11] . 
He uses the relaxation time approximation and finds that C^i n (t) = i~ 3 / 5 for large 
t. Using a resolvent expansion, in [13] we prove corresponding sharp bounds and 
thereby confirm the relaxation time approximation in this particular case. For a 
finite chain of length N with thermal reservoirs at both ends, the energy transport 
is then anomalous and the thermal conductivity diverges as iV 2//5 , which seems to 
be in agreement with molecular dynamics. For a more detailed discussion we refer 
to [TH], Section 6. 

5 Gaussian fluctuation theory 

Energy transport can be viewed in the more general context of time-dependent 
Gaussian fluctuation theory close to thermal equilibrium. For low density gases this 
link is reviewed in [16] with further examples discussed in [17]. The purpose of this 
section is to explain how phonon kinetic theory makes no exception. In [TH1 EZ] 
spatial variation is included. Since our exposition deals only with the spatially 
homogeneous system, we stick to such a set-up also for the fluctuation theory. 

Physically, one considers time-dependent fluctuations in equilibrium for the num- 
ber of phonons with wave number k. Technically one has to sum over phonons in 
a small volume element in fc-space. To be more precise we partition the tours 
T = [—1/2, 1/2] by a grid with spacing e and denote it by T e . (T £ ) d corresponds to 
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the crystal volume [1, . . . , l] d C 1> d with periodic boundary conditions, I — 1/e. Let 
/ : T d — > K be a smooth test function. Then the fluctuation field, indexed by / and 
t, is defined through 

W, t) = e d ' 2 Y, f( k ) *) - (a £ (kya%k))p) . (5.1) 

k£(T e ) d 

a £ (k,t) depends on e through the finite crystal volume e~ d , through setting A 2 = e, 
and through the rescaled time e~ l t in microscopic units. The claim is that, in 
distribution, the limit 

lim£ e (/, *)=&(/) (5.2) 

exists and that the limit random field is classical. In fact, the limit field should 
be jointly Gaussian and governed by the linear Langevin equation 

^ t (k) = A^(k) + B Vt (k) , (5.3) 

where £*(/) = jfddkf{k)E st {k). A is the generator from the linearized Boltzmann 
equation, compare with f)3.20p . and rj t is normalized Gaussian white noise with 

E( Vt (k) Vt/ (k')) = 5{t - t')5(k - k') . (5.4) 

The linear operator B controls the strength and correlations for the noise input to 
the various fc-modes. 

The main observation of the fluctuation theory is the relationship between A and 
B through the equal-time equilibrium fluctuations. We set, as a linear operator, 

(g, Cf) = lim(r(^ 0)Of, 0))p. (5-5) 

e— >0 

Using that 

(a*(k 1 )a(k 2 )a*(k 3 )a(k4))p - (a*{k 1 )a(k 2 ))p(a(k 3 )*a(k4))p 

= 5(h - k 4 )5(k 2 - h)Wp(h)W p (k 2 ) , (5.6) 

one obtains 



(g,Cf)= dkg(k)W p (k)Wp(k)f(k), (5.7) 

in other words C is the operator of multiplication by WpWp. The fluctuation- 
dissipation relation takes then the form 

AC + CA* = ~BB* . (5.8) 

Since AC = L = L* = CA*, one concludes that the noise strength is 

BB* = 2L. (5.9) 

A posteriori this identity explains also the at first sight unexpected linearization 
in (I3.18p . Only then the linearized operator is symmetric, as is obvious from (15.91) . 
Solving (15. 3p . the covariance of the stationary fluctuation field is given by 

Ug)Uf)) = (g, e Altl cf), (5.io) 

in agreement with the special case / = g = (2tt)~ 1 (£- Vuj)uj of interest in Section[3l 
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